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Abstract 
We present the first determination of the visible reflection spectrum 
of a dark spot in Neptune’s atmosphere, ‘NDS-2018’, analogous to the 
‘Great Dark Spot’, discovered by Voyager 2 in 1989. Observations were 
made in 2019 with the MUSE Integral Field Unit spectrometer at the 
Very Large Telescope. We show that the feature is caused by a dark- 
ening at wavelengths shorter than 700 nm of a deep layer of aerosol 
in Neptune’s atmosphere at a pressure of ~5 bar, which we suggest 
is coincident with the main H2S condensation layer. In addition, we 
have detected an apparently short-lived companion deep bright spot 
feature, ‘DBS-2019’, on the south-west edge of NDS-2018, with a spec- 
tral signature consistent with a brightening of the same 5-bar layer at 
wavelengths longer than 700 nm. This feature is fundamentally different 
from previously studied companion clouds of the Voyager-2 Great Dark 
Spot, which were located much higher in the atmosphere at 0.6-0.2 bar. 


1 Introduction 


Planetary-scale vortices are commonplace features of giant planet atmospheres. 
The most famous example is Jupiter’s Great Red Spot (GRS), whose red colour 
is explained by the presence of a blue-absorbing ‘chromophore’/1] in a haze 
above an anticyclonic vortex with a cold core above the mid-plane at ~2-5 
bar and a warm core below[2, 3]. The GRS lies at a latitude of ~22°S and is 
currently ~15,000 km wide. While the GRS has been observed since at least 
the middle of the 19th Century, no comparably large and visible vortex had 
been seen in another giant planet atmosphere until the arrival of Voyager 2 at 
Neptune in 1989/4]. During the approach of this spacecraft, a large, dark anti- 
cyclonic vortex, the ‘Great Dark Spot (GDS)’, was seen by the Imaging Science 
Subsystem (ISS) at ~20°S with a size of ~10,000 km. This spot was dark at 
blue wavelengths, but became indiscernible at wavelengths longer than 700 
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nm. Although it is tempting to think this feature was similar to Jupiter’s GRS, 
it appeared to be mostly located at deeper levels, although it was associated 
with ‘companion’ methane ice condensation clouds in the upper atmosphere 
(~0.6-0.2 bar). The GDS was observed for about 7 months by Voyager 2, but 
was never seen again and thus appears to have been a short-lived feature[5]. 
Furthermore, since Voyager did not have the capability to provide more than 
filter-imaging observation of the reflected-sunlight spectra of this vortex its ver- 
tical structure remained largely unknown. However, since the discovery of the 
GDS, several more short-lived dark spots in Neptune’s atmosphere have been 
detected in Hubble Space Telescope (HST) filter-imaging observations with 
the Wide Field and Planetary Camera 2 and Wide Field Camera 3 (WFC3) 
instruments[6—8], in both the northern and southern hemispheres. The most 
recent example is a northern hemisphere dark spot, discovered in 2018 at 23°N 
(NDS-2018)[9]. This spot had a similar size to the GDS and was then seen to 
drift equatorwards[10], before apparently disappearing in late 2022[11]. Nep- 
tunian dark spots are characterised by low reflectance at short wavelengths 
(A < 700 nm), but have not been detected at longer wavelengths [4, 7]. 

A recent reanalysis of multiple space- and ground-based observations from 
0.3 — 2.5 um[12], hereafter ‘IRW22’, found that the spectra of both Uranus 
and Neptune can be modelled very accurately with a simple atmospheric 
aerosol structure comprised of three main layers: 1) a deep H2S/photochemical- 
haze aerosol layer with a base pressure > 5-7 bar (Aerosol-1); 2) a layer of 
methane/photochemical-haze just above the methane condensation level at 
1-2 bar (Aerosol-2); and 3) an extended layer of small photochemical haze 
particles extending into the stratosphere (Aerosol-3). IRW22 also analysed 
HST/WFC3 observations of NDS-2018[9], and Voyager-2/ISS observations of 
the GDS[4] together with a dark band near 60°S, dubbed the ‘South Polar 
Wave’ (SPW)[13]. IRW22 deduced that dark features are most likely caused 
by a darkening of the deep Aerosol-1 layer, although a ‘clearing’ of this layer 
could not be ruled out. The deep location (p > 3 bar) of the SPW was also 
deduced from a previous analysis of HST filter-imaging observations from 1994 
to 2008/14]. Unfortunately, since all previous observations of Neptune’s dark 
features have been restricted to a few filter-imaging wavelengths it has not 
been possible to say definitively whether these features are dark due to the 
presence of chromophores, like Jupiter’s GRS, whether they are due to cloud 
‘clearing’, or whether they are darkened by a completely different mechanism. 
Hence, the goal of our study was to record a complete visible/near-IR. spec- 
trum of a dark feature and determine if the spectral ‘signature’ obtained could 
be used to differentiate between these different darkening scenarios. 


2 Results 


Observations and Processing. As part of a global collaboration to observe 
and understand the NDS-2018 feature[10], we observed Neptune on 17th/18th 
October and 13th November 2019 with the Multi Unit Spectroscopic Explorer 
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(MUSE) Integral Field Unit (IFU) spectrometer at the Very Large Telescope 
(VLT) of the European Southern Observatory at La Paranal, Chile. The MUSE 
hyperspectral instrument records images at 2048 wavelengths simultaneously 
from 475 to 993 nm over a 7.5” x 7.5” field of view in Narrow-Field Mode. 
These data were smoothed to resolution of 2 nm (See Methods) to improve 
the SNR, but even with the GALACSI[15] adaptive optics system active NDS- 
2018 was barely distinguishable in the raw data, requiring the development of a 
tailored deconvolution technique (see Deconvolution Methods in the Methods 
Section). Figure la shows the observed appearance of Neptune in our clear- 
est, smoothed MUSE dataset, ‘Obs-6’, including NDS-2018 (Table 1) at three 
different wavelengths: 1) 551 nm, which sounds deep in the atmosphere and is 
most sensitive to dark features; 2) 831 nm, which also sounds deep in the atmo- 
sphere, but is not strongly affected by Rayleigh scattering; and 3) 848 nm, at 
the edge of a strong methane absorption band, which is most sensitive to reflec- 
tion from haze high in the atmosphere (p < 0.5 bar). These wavelengths were 
chosen for being spectrally representative and relatively free from artefacts (see 
Methods). The three columns in Fig. la are: 1) non-deconvolved appearance; 2) 
deconvolved appearance; and 3) contrast-stretched and flattened deconvolved 
appearance (i.e., corrected for disc-averaged limb-darkening). Although invis- 
ible in the non-deconvolved data, NDS-2018 is visible at the top right of the 
551-nm deconvolved images (at 15°N, 10°E of the central meridian), but is 
invisible at the longer wavelengths shown here. Similarly, the dark SPW[13] at 
~60°S is visible at 551 nm, but not at longer wavelengths. At 831 nm, where 
we can detect reflection from cloud/haze lying very deep in Neptune’s atmo- 
sphere (5—7 bar), we see bright zones and dark belts with a width of ~ 20° and 
a particularly bright zone at a latitude of 65—70°S, just south of the SPW and 
just north of the frequently observed South Polar Features (SPF)[4, 16]. More 
intriguingly, a bright spot appears just to the south west of NDS-2018 at 10°N, 
0°E, which we refer to as ‘Deep Bright Spot 2019’ (DBS-2019). We found the 
NDS-2018 and DBS-2019 features in all five dark spot observations (Table 1) 
and saw both objects rotating eastwards with the planetary rotation, proving 
these are real features and not systematic artefacts (Supplementary Fig. 4). No 
new features are visible in the 848-nm image, although the bright clouds near 
the left edge of the disc, which we shall refer to as Shallow Bright Spots (SBS), 
spanning ~70°S to ~30°S, are more clearly defined; the spot near ~70°S could 
be more properly described as a South Polar Feature (SPF), first identified by 
Voyager 2[4]. The spatial relationship between the NDS-2018 and DBS-2019 
features can be seen more clearly Fig. 1b, which shows a false-colour compos- 
ite of the deconvolved and enhanced appearances, with DBS-2019, appearing 
red, lying just to the south west of the darker NDS-2018. 

Fig. 1c shows the difference between the spectra of the discrete features and 
those expected at their locations from a latitudinal-average Minnaert limb- 
darkening analysis (see Methods). The NDS-2018 difference spectrum shows 
increasing darkening relative to the expected background at wavelengths < 700 
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Table 1 VLT/MUSE Neptune observations (Narrow-Field Mode) containing the 
NDS-2018 feature . 


Obs. ID. Date Time (UT) exposure time airmass 
6 October 18th 2019 00:01:20 120s 1.218 
7 October 18th 2019 00:18:08 120s 1.172 
8 October 18th 2019 00:22:53 120s 1.161 
9 October 18th 2019 00:27:36 120s 1.158 
10 October 18th 2019 00:32:19 120s 1.141 


N.B., the Observation ID simply relates to order of observation in our data set. 


nm, with strong absorption bands of gaseous methane absorption visible, indi- 
cating the darkening is located at considerable depth in Neptune’s atmosphere. 
In contrast, the DBS-2019 feature has minimal signature at short wavelengths, 
but has well defined, narrow increases in reflectance near 750, 830 and 930 nm. 
These wavelengths coincide with methane absorption ‘windows’ and the nar- 
row reflectance peaks strongly indicates DBS-2019 also to be a deep feature, 
formed by reflectance changes at pressures > 3 bar (Supplementary Figs. 8 
and 9). The difference spectra of the Shallow Bright Spots (SBS) are very dif- 
ferent from DBS-2019 and are visible at all wavelengths longer than 600 nm, 
indicating these to be high (0.2-0.6 bar) methane ice clouds. No other deep 
bright spots like DBS-2019 were seen in any of our other observations at any 
latitude. 

Radiative Transfer Analysis of NDS-2018 and DBS-2019. We 
modelled the spectral signatures of NDS-2018 and DBS-2019 using the 
NEMESIS[17-19] radiative transfer model, adapting the procedure used by 
IRW22[12] for the analysis of HST/STIS, IRTF/SpeX and Gemini/NIFS obser- 
vations of Uranus and Neptune from 0.3 to 2.5 um. Updating the procedure 
slightly from IRW22 we parameterised the lowest ‘Aerosol-1’ haze as a ver- 
tically thin layer, centred at 5 bar, but left the other model parameters 
unchanged (see Methods). We fitted both the observed spectra at the discrete 
cloud feature locations, and also the expected spectra at the same locations 
reconstructed from latitudinally-averaged Minnaert fits (see Methods). As 
IRW22 found the limb-darkening of the reflectivity of the Voyager-2 Great 
Dark Spot and South Polar Wave cannot be explained by changes in opacity 
or reflectivity of the Aerosol-2 or Aerosol-3 layers, and since the shape of the 
DBS-2019 difference spectrum also suggests perturbations at depth, we limited 
our analysis of both features to perturbations of the Aerosol-1 layer only. The 
measured spectrum of the features and the expected spectrum at the same loca- 
tion determined from our latitudinally-averaged limb-darkening analysis were 
thus compared to determine any necessary corrections to the Aerosol-1 opac- 
ity and scattering properties, keeping all other atmospheric properties fixed to 
the Minnaert-analysis of the respective latitudinal band (See Methods). 

The NDS-2018 difference spectrum and our fits to it are shown in Fig. 2a 
(the individual NDS-2018 and expected background spectra and their respec- 
tive fits are shown in Supplementary Fig. 5). The NDS-2018 dark spot is 
darker than the expected background at wavelengths shorter than 700 nm, 
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but has little signature at longer wavelengths. We found that simply reduc- 
ing the Aerosol-1 opacity from that found for the background provided a poor 
fit to the difference spectrum. However, we could achieve a very good fit to 
the difference spectrum when, in addition to slightly lowering the Aerosol-1 
opacity, we also allowed the imaginary refractive index spectrum (Nimag) of 
these particles to vary (Fig. 2b), increasing Nimag (and thus reducing their 
single-scattering albedo, making them less reflective) in the 500 — 600 nm 
range. Figure 2b compares the fitted imaginary refractive index spectrum of 
the Aerosol-1 particles with those expected for all other aerosol types at the 
NDS-2018 location, while Fig. 2c shows the fitted aerosol profiles at 10-20°N 
latitude band and the inferred perturbations to the Aerosol-1 profile. For ref- 
erence, Supplementary Fig. 7 compares the imaginary refractive index spectra 
shown in Fig. 2b with those deduced for the blue-absorbing chromophore in 
Jupiter’s atmosphere[20, 21], showing the Aerosol-1 particles to have a very 
different spectral shape to those seen in Jupiter’s atmosphere. Possible candi- 
dates for the Neptune ‘chomophore’ are discussed by IRW22[12]. We conclude 
that the complete spectral signature of NDS-2018 captured by VLT/MUSE 
allows us to rule out the cloud ‘clearing’ scenario for dark spots with high 
confidence and demonstrates that dark spots are mostly the result of changes 
in the properties of the Aerosol-1 particles that reduce their single-scattering 
albedo and make them less reflective at short wavelengths. 

Turning to DBS-2019, the bright, narrow difference peaks at 750, 830 and 
930 nm (Fig. 3a) suggest a thickening or brightening of the aerosols in the deep 
atmosphere and Supplementary Figs. 8 and 9 show the differences to the com- 
puted spectra when an additional scattering layer is introduced at different 
pressure levels. We find that an additional haze layer at pressure <~4 bar pro- 
duces difference features at wavelengths longer than 700 nm that are too wide, 
while an additional layer at pressure >~6 bar would need to be enormously 
thick in order to produce an observable signature. Hence, the DBS-2019 and 
expected background spectra were again compared to determine any neces- 
sary corrections to the opacity and scattering properties of the Aerosol-1 layer, 
keeping all other atmospheric properties fixed to the Minnaert-analysis of the 
5 — 15°N latitudinal band. Our fitted difference spectra are shown in Fig. 3a, 
while fits to the individual spectra are shown in Supplementary Fig. 6. Again, 
we found that we could fit the difference spectrum only by modifying both 
the opacity and the optical properties of the Aerosol-1 particles, here finding 
it necessary to reduce the imaginary refractive index at longer wavelengths to 
increase their single-scattering albedo and make them more reflective. 

Although our variable njmag model fits the DBS-2019 bright spot very 
well at longer wavelengths, it was slightly too reflective at shorter wavelengths 
(~ 2%). One way we found to correct this was to increase also the mean 
radius of the Aerosol-1 particles from 0.1 to 0.7 um. However, this correction 
is not unique and depends on our assumed particle size distribution, not just 
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the mean radius of the aerosols. In addition, the close proximity of DBS- 
2019 to NDS-2018 means we cannot exclude the possibility of some cross- 
contamination of the spectra. We conclude that the DBS-2019 feature is caused 
by a brightening of the Aerosol-1 particles at ~5 bar at wavelengths longer 
than 700 nm, which is possibly combined with a change in the particle size 
distribution. 


3 Discussion 


We find that the darkness of NDS-2018 is caused by a chromophore that makes 
the particles in a deep Aerosol-1 layer (IRW22[12]) at ~5 bar less reflective at 
wavelengths shorter than 700 nm. Here we consider two possible explanations 
for this darkening. 

Firstly, it is possible that a dark chromophore is being introduced to the 
Aerosol-1 layer, either from above or below. We think the former scenario is 
unlikely since we have good sensitivity to aerosols at pressures less than 5 
bar and would easily detect such particles as they fell down to the Aerosol-1 
level. The latter scenario is harder to rule out since we cannot see to deeper 
levels through the Aerosol-1 layer. Hence, dark material could be upwelling 
from below the dark spot, although what process might produce such dark 
material is unknown. One hypothesis is that increased abundances of a light- 
sensitive gas, e.g., H2S, might be upwelling and photolysed by ultra-violet 
(UV) radiation, leading to a dark, photochemical product. From the best-fit 
Neptune solution of IRW22[12], we find that although the overlying Aerosol- 
3 and Aerosol-2 aerosols are predicted to absorb some of the UV flux, in the 
absence of a significant opacity of Aerosol-1 particles considerable UV could 
scatter conservatively down to very deep levels and photolyse light-sensitive 
gases (Supplementary Figs. 10 and 11). The products of this photolysis might 
then increase the UV absorption of the Aerosol-1 particles, thus regulating the 
photolysis and also increasing the UV/visible darkness of the particles. The 
vertical motion of air inside dark spots is unknown, but will likely depend 
on the mid-level of the vortex, which could potentially be at the methane 
condensation level for Neptune[12], leading to downwelling at 5 bar. More 
significantly there may be very little mixing between the air inside and outside 
of the anticyclonic vortex and so particles within the vortex may have time to 
be considerably more photo-processed (and so darkened) than those outside. 

A second scenario for dark spot coloration is that the Aerosol-1 particles 
are formed from H2S ice condensed on to darker, photochemically-produced 
haze particles transported down from the stratosphere (IRW22). IRW22 sug- 
gest that dark spots may be caused by local heating at the ~5-bar level in the 
deep half of an anticyclonic vortex, subliming the H25 ice to reveal their darker 
photochemical haze cores. Such a model of dark spot formation would lead to a 
reduction in the Aerosol-1 particle sizes and we do find that the Aerosol-1 par- 
ticles in the dark spot have both lower scattering albedo at short wavelengths 
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and also slightly lower opacity, which is consistent with this hypothesis. Sub- 
limation of HS might also lead to the formation of more chromophore if the 
UV flux is sufficiently high, which would darken the particles even more. How- 
ever, H2S sublimating at 5 bar could re-condense at lower pressure if the air is 
upwelling, but we see no evidence of any such recondensation in the retrieved 
aerosol structure. This could mean that the air is indeed downwelling but it 
could also mean that the hypothesis is incorrect. Clearly, more work is needed 
to determine whether this scenario is plausible. Unfortunately, such an analy- 
sis is hampered both by our lack of knowledge of the complex refractive indices 
of either constituent and also by the lack of any measurements to test whether 
the proposed thermal increases within the dark spot are realistic. 

Neither scenario for the nature of dark spots and what causes their 
coloration bears much resemblance to Jupiter’s Great Red Spot, another 
anticyclonic vortex, which extends right up to the tropopause, is reddened 
by a blue-absorbing chromophore in its upper aerosol layers and shows ele- 
vated abundances of aerosols and phosphine[22], and a slight elevation in 
ammonia|22, 23]. Instead, Neptune’s dark spots appear to be deep and do not 
show elevated abundances of tropospheric methane, which would have been 
detected in our MUSE observations. Work is ongoing to determine if either 
the stratospheric methane abundance or tropospheric hydrogen sulphide abun- 
dance is perturbed over NDS-2018, through analysis of VLT/SINFONI and 
Keck/OSIRIS observations. 

The fact that the NDS-2018 lies so close to the adjacent DBS-2019 fea- 
ture, which seems to be caused by a spectrally-dependent brightening of the 
aerosols in the very same Aerosol-1 layer at ~5 bar, is intriguing and suggests 
some sort of connection. Voyager 2 images showed that the dark oval, DS2 
(at 55°S), developed bright clouds at its centre that changed on timescales 
of hours. Their morphology, rapid variability and brightness suggested a con- 
vective origin, distinct from the companion clouds of the GDS/4, 16]. It could 
be that DBS-2019 is a small region of localised upwelling at the edge of the 
vortex, perhaps increasing H2S condensation there and making the particles 
more reflective at longer wavelengths. Whatever DBS-2019 is, it seems to be 
an ephemeral feature. Although we need the high spectral resolution of MUSE 
to differentiate between deep features such as DBS-2019 and the more com- 
mon and higher (~0.6-0.2 bar) methane ice cloud Shallow Bright Spots (SBS), 
we find that the DBS-2019 feature would have been visible in the F845M filter 
had HST/WFC3 been observing at the same time as VLT/MUSE. However, 
although HST/WFC3 observations of Neptune were made on 28/29th Septem- 
ber 2019, less than three weeks earlier, no trace of the DBS-2019 feature was 
seen. Furthermore, less than a week before our observations, Keck/OSIRIS 
(H-band IFU) detected neither the NDS-2018 nor the DBS-2019 on 12/13th 
October. This non-detection of the dark spot at H-band wavelengths was 
not unexpected, but since DBS-2019 is interpreted here as being caused by 
a spectrally-dependent brightening and possibly thickening of the Aerosol-1 
haze, it is interesting that it was not seen in the Keck/OSIRIS data. This could 
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mean either that DBS-2019 was not present just a week before we observed 
with VLT/MUSE, or alternatively that its signature is insignificant at longer 
wavelengths. 

Other ephemeral, short-lived spots have been seen before in Neptune’s 
atmosphere, such as the South Polar Feature (SPF), a group of small, tran- 
sient clouds originally observed from 68°S to 75°S by Voyager 2[16, 24]. Such 
features have since been seen in H-band observations by Keck/NIRSPEC[25], 
Keck/NIRC2[26] and Gemini/NIFS[27], and may be linked with the South 
Polar Wave (SPW)[14]. Most SPF clouds seem to reside in the upper tropo- 
sphere, just below the tropopause at ~0.2 bar, which suggests they are usually 
similar to our Shallow Bright Spots (SBS). However, Gemini/NIFS observa- 
tions made in 2009/27], and perhaps Keck/NIRC2 images from 2003[26], found 
some SPF clouds to lie deep in the atmosphere at p > 1 bar, and to have a life- 
time of only 1-2 days, similar to the DBS-2019 feature. Since this latitude lies 
next to the dark South Polar Wave (SPW) at ~60°S, there would again appear 
to be a possible link between deep, bright clouds and dark, visible-wavelength 
features. 

The NDS-2018 dark spot studied here appears to be very different from 
the Great Red Spot anticyclonic vortex seen on Jupiter, and is caused by a 
spectrally-dependent darkening of the particles in a layer at ~5 bar. This shows 
that the appearance of giant planet anticyclones depends critically on their 
background environment — i.e., where the midplane of the anticyclone sits with 
respect to vertical features such as the tropopause and condensation cloud 
levels. To further probe these differences requires thermal measurements to be 
made at multiple depths in the atmosphere, ideally from a future proposed 
space mission to either Neptune or Uranus, combined with better laboratory 
determinations of the refractive index spectra of potential aerosols and more 
detailed numerical dynamical simulations. 


4 Methods 


Observations. MUSE is an Integral Field Unit Spectrometer, where each 
pixel of its 300 x 300 pixel Field of View (FOV) is a complete spectrum 
covering 2048 wavelengths from 475 to 933 nm at a spectral resolution of 2000 
— 4000, thus forming a ‘cube’ of data. MUSE was operated in its Narrow- 
Field Mode, which has a projected FOV of 7.5” x 7.5” and an individual 
‘spaxel’ size of 0.025” x 0.025”. Since the spectral resolution of our analysis was 
limited to the available sources of methane absorption data[28], we spectrally 
averaged the MUSE cubes with a triangular instrument line shape of Full- 
Width-Half-Maximum (FWHM) 2 nm, and sampled every 1 nm to achieve 
Nyquist sampling. This reduced the number of wavelengths in the MUSE cube 
from 2048 to 459 and also improved the signal-to-noise ratio. 

Observations were made using the GALACSI adaptive optics system[15] 
with a laser guide star, which is quoted to achieve a spatial resolution of 0.05 
— 0.08”, much smaller than the expected horizontal size of NDS-2018 of ~0.2”, 
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based on the HST/WFC3 observations. We used HST observations[9, 10] to 
estimate when the feature would be visible in our data and found that five 
observations were made when NDS-2018 was probably close to the central 
meridian on 18th October 2019: these are listed in Table 1. Of these five obser- 
vations, the one with the best spatial resolution was Observation 6 (‘Obs-6’), 
for which NDS-2018 at 15°N, was predicted to be 3.3°E of the central merid- 
ian. While the GALACSI adaptive optics system achieves a spatial resolution 
of ~0.06” at 800 nm in Narrow-Field Mode, this performance deteriorates to 
~0.2” at shorter wavelengths (Supplementary Fig. 1). Since the dark spot con- 
trast increases at shorter wavelengths [4, 6, 7], but spatial resolution decreases, 
we found that wavelengths near 551 nm were optimal for dark spot detec- 
tion. However, even at these wavelengths, we were unable to clearly distinguish 
the small, faint NDS-2018 from its surroundings in the raw data (Fig. la). 
Hence, we developed a deconvolution algorithm to improve the spatial res- 
olution still further, using an approach based on a modified version of the 
CLEAN algorithm[29]. Here, instead of operating on the single brightest pixel 
of an image in each iteration, as in the original implementation of CLEAN, 
our scheme operates on all pixels above some fractional threshold value[30] 
(See Deconvolution Methods later in this section). The wavelength-dependent 
point spread function (PSF) used during this deconvolution was obtained from 
a standard star observation made 20 minutes before the Neptune observation 
(Supplementary Table 1), during very similar atmospheric conditions, and the 
form of this PSF at 551 nm is shown in Supplementary Fig. 2. An example 
test deconvolution of a synthetically generated and convolved image is shown 
in Supplementary Fig. 3, demonstrating the efficacy of this technique. This 
deconvolution algorithm was applied to all wavelengths of our smoothed cube 
and leads to a dramatic improvement in spatial resolution. 

In Fig. 1a we show individual, smoothed wavelengths, rather than an aver- 
age of several nearby to indicate the data quality, and we picked representative 
wavelengths to show the key features. For the first row of Fig. la, we could 
easily have shown images at a neighbouring wavelength (e.g., 550nm or 552 
nm) rather than 551 nm, but 551 nm was chosen as being representative of 
that wavelength region, showed the dark spot clearly, and was relatively free 
from artefacts. Similarly, for the second row of Panel B, we chose the 831-nm 
image as it is near the centre of one of the DBS-2019 reflectivity peaks (Fig. 
1c) and shows this feature clearly. Finally, for the methane-absorbing image 
(third row of Fig. 1a) we chose the 848-nm image since this is less noisy than 
observations in band centre and is sensitive to slightly deeper levels, revealing 
the haze band near the equator. 


Minnaert Limb-darkening Approximation. The reflectance spectra of 
Neptune are found to be well approximated by the Minnaert limb-darkening 
approximation[31] in the MUSE spectral range[32]. Here, the reflectivity, 
I, of an observation at a particular wavelength may be approximated as 
I = Ippkuk~, where pu and po are the cosines of the viewing and solar zenith 
angles, respectively, and Jp is the fitted nadir reflectance and k is the fitted 
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limb-darkening parameter. For each wavelength in the wavelength-smoothed 
MUSE data, the reflectivities across the disc were analysed in latitudinal bands 
of width 10° (spaced every 5° to achieve Nyquist sampling), and the fitted 
Minnaert parameters, Io(¢, A) and k(¢, A) recorded, where ¢ is the mean lati- 
tude and A is the wavelength. These coefficients were then used to generate the 
expected spectra at two observing zenith angles for each latitude band that 
could be fitted simultaneously with our radiative transfer model NEMESIS[17], 
following IRW22[12]. 


Radiative Transfer Modelling. We modelled these observations using 
the NEMESIS[17-19] radiative transfer and retrieval model, following a pre- 
vious analysis of HST/STIS, IRTF/SpeX and Gemini/NIFS Uranus and 
Neptune observations from 0.3 to 2.5 um, IRW22[12]. As IRW22[12] note, these 
retrievals are difficult and prone to degenerate solutions since we have very 
little a priori knowledge of the vertical structure of clouds, the atmospheric 
composition, or the scattering properties of the cloud and haze particles. 
IRW22 reduce this degeneracy by analysing a wide wavelength range and sev- 
eral viewing angles simultaneously and the reader is referred to this paper 
for a full discussion. IRW22 found that the spectra of both planets were well 
modelled with three main aerosol layers: 1) ‘Aerosol-1’, a deep aerosol layer 
with a base pressure > 5-7 bar, assumed to be composed of a mixture of H2S 
ice and photochemical haze, of assumed mean radius 0.05 um; 2) ‘Aerosol-2’, 
a layer of photochemical haze and methane ice, of mean radius ~ 0.5 um, 
confined within in a layer of high static stability at the methane condensation 
level at 1-2 bar; and 3) ‘Aerosol-3’, an extended layer of small photochemical 
haze particles (r ~ 0.05 um), probably of the same composition as the 1—-2-bar 
layer, extending from this level up through to the stratosphere. For Neptune, 
an additional thin layer of micron-sized methane ice particles at ~0.2 bar 
was required to explain enhanced reflection at methane-absorbing wavelengths 
longer than 1.0 wm, which was not required in this analysis. 

We initially fixed the complex refractive index spectra to the best-fitting 
solutions found from the IRW22[12] analysis of their combined 0.3-2.5 um 
Neptune observations, and similarly fixed the mean radius of the Aerosol- 
2 particles to the average best-fitting value of 0.7 wm. We fixed the mean 
radius of the Aerosol-1 particles to 0.1 um (increased from the IRW22 estimate 
of 0.05 um to be more in line with expectations for condensed fog parti- 
cles), but kept the mean radius of the Aerosol-3 particles fixed at 0.05 ym. 
Then, for each 10°-wide latitude bin (stepped every 5° in latitude to achieve 
Nyquist sampling) we fitted to spectra reconstructed from the Minnaert fit- 
ted parameters for pure back-scattering at solar and viewing zenith angles 
of 0° and 61.45°, respectively. Two angles were chosen like this so that we 
could simultaneously fit to both the mean reflectivity and the observed limb- 
darkening/limb-brightening, and these two particular angles were chosen as 
they coincided with two of the zenith angles used in the zenith angle quadra- 
ture scheme of our plane-parallel multiple scattering model|[33]. The random 
errors on these reconstructed spectra were smaller than the modelling error of 
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our forward model and so, following IRW22[12] the uncertainties on the spec- 
tra were set to 1/50 of the maximum nadir reflectivity within 50 nm of each 
wavelength[12]. Following IRW22, these errors were further reduced by a fac- 
tor of 2 between 800 and 900 nm to force a closer fit in this region, which is 
most sensitive to the methane mole fraction. 

We used the whole 475 — 933 nm range of the MUSE data, except for 
the 578 — 605 nm range reserved for the laser guide star, and sampled at 
the same set of wavelengths as previously used by the IRW22[12] analysis 
of the HST/STIS observations. We also included additional wavelengths at 
wavelengths lower than those observed by MUSE to ensure that the Raman- 
scattering[34] components of these shorter wavelengths were approximately 
accounted for in our MUSE simulations. This additional part of the spectrum 
was reconstructed from HST/STIS observations that were scaled to match 
the MUSE observations at 475 nm. The uncertainties on these added/scaled 
HST/STIS spectral points were set to 100% so that our fitting model would 
not attempt to fit them closely, but would still include their Raman-scattered 
contribution to the MUSE wavelengths. 

Although this standard reference model fitted the VLT/MUSE data well 
in general, there were several deficiencies. Firstly, the deconvolved MUSE data 
have rather higher spatial resolution than the HST/STIS data and thus there 
is more limb-brightening at methane-absorbing wavelengths than is seen in 
the HST/STIS data. Increased limb-brightening implies haze particles that 
are more scattering than those found to be consistent with HST/STIS. Hence, 
we Minnaert-analysed the entire disc of Neptune (masking out discrete cloud 
features), reconstructed disc-averaged spectra at 0° and 61.45° zenith angles 
and re-retrieved the global-mean cloud opacities, methane mole fraction and 
the imaginary refractive index (Nimag) spectra of the three aerosol types. Sec- 
ondly, while the reference model fitted the disc-averaged data well, it was 
less successful in fitting the latitude-dependence of these MUSE data, with 
cross-correlation observed between the retrieved Aerosol-1 opacity, Aerosol-1 
fractional scale height, Aerosol-2 opacity and methane abundance. To reduce 
this degeneracy we simplified the Aerosol-1 parameterisation to be a single 
thin layer fixed at 5 bar with no overlap with Aerosol-2, rather than using 
the previous model with a fixed base pressure and variable scale height that 
led to both Aerosol-1 and Aerosol-2 particles at some altitudes and some 
degeneracy. We found this revised model matched the data equally well and 
was less degenerate. We will return to latitudinal changes in a forthcoming 
paper currently in preparation. Making this change to the Aerosol-1 profile 
parameterisation we further updated the retrieved global-mean atmospheric 
properties and Nimag Spectra of the three Aerosol types from the disc-averaged 
Minnaert-reconstructed spectra at 0° and 61.45° zenith angle. 

To analyse the NDS-2018 and DBS-2019 spectra, we concentrated on the 
5 — 15°N and 10 — 20°N bands respectively. In each band, initially fixing the 
Nimag Spectra of the three aerosol types to the values determined from the 
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disc-averaged Minnaert analysis, we fitted the latitudinally-averaged Minnaert- 
reconstructed spectra at 0° and 61.45° zenith angle to determine the opacities 
of the three aerosol components, the pressure of the Aerosol-2 layer and the 
methane deep abundance. The spectra of the NDS-2018 and DBS-2019 were 
averaged over a 2x 2 box of pixels centred on the features. To properly compare 
these spectra with those expected at these locations, the Minnaert-parameters 
in the respective latitude bands were used to reconstruct the spectra that we 
would to see at the NDS-2018 and DBS-2019 locations, which we call the 
‘background’ spectra. To fit these background spectra as closely as possible, we 
found it necessary to refit the Aerosol-1 opacity and Aerosol-1 nimag spectra, 
keeping all other parameters fixed at their latitude-averaged values. The NDS- 
2018 and DBS-2019 spectra were then fitted in the same way. In all cases the 
uncertainties on the spectra were set to the same forward-modelling values as 
previously described. Supplementary Figs. 5 and 6 show our best fits to the 
NDS-2018 and DBS-2019 spectra, and their respective background spectra, 
and the differences between these fits are the difference spectra shown in Figs. 
2 and 3. The error ranges on the measured difference spectra shown in Figs. 2 
to 3, set to the standard deviation over the respective 2 x 2-averaging boxes, 
are shown for reference and were not used in the fitting process. 


Scattering Properties. We retrieved the imaginary refractive index spec- 
tra Nimag for all the aerosols and then reconstructed the nreq spectra using a 
Kramers-Kronig analysis, assuming that all had the same real refractive index 
of 1.4 at a reference wavelength of 800 nm. The nimag spectra were tabulated 
every 100nm, and a correlation length of 100 nm assumed to achieve some 
degree of wavelength smoothing. From the derived complex refractive index 
spectrum we then calculated the extinction cross-section and single scatter- 
ing albedo spectra using Mie theory. However, since the particles in the Ice 
Giant atmospheres will be ices, rather than liquid, we approximated the Mie- 
calculated phase functions with combined Henyey-Greenstein phase functions, 
which average over features peculiar to spherical particles, such as the ‘glory’ 
(i.e., the increased scattering seen exactly 180° away from the incident beam at 
the antisolar point), and the coloured ‘rainbow’, commonly seen 40-50° away 
from the antisolar point. The derived particle scattering properties used in our 
fits to the NDS-2018, DBS-2019 spectra, and their respective backgrounds, are 
listed in Supplementary Tables 2 — 7. 


Deconvolution Methods. We used a modified version of the CLEAN 
algorithm[29], which is frequently used in the radio telescope community [35], 
but not often in the optical [36]. Similarly to [30], instead of single-point sub- 
tractions, pixels above some selection threshold, ts, were convolved with the 
instrumental PSF and a fraction of the result, the loop-gain (gioop), subtracted 
from the “dirty map” each iteration: we will refer to this as MODIFIED-CLEAN. 
We used the standard-star observations taken as part of the observing run as 
our PSF, see Supplementary Table 1 for details. In our case ts was determined 
dynamically by choosing an initial threshold, ti, from successive applications 


Springer Nature 2021 TeX template 


14 Cloud Structure of Dark Spots and Storms in Neptune’s Atmosphere 


of ‘Otsu’ thresholding [37], and applying a user-supplied factor, tractor (set to 
0.3 for this work), to the interval between t; and the data maximum such that 


ts = (data maximum) x tractor + (1 — tractor) X ti. 


Therefore, tractor is equivalent to the trim contour described in [30], and our 
choice of t; from many candidate Otsu thresholds determine which features 
are deconvolved first. 

Otsu thresholding [38] is a method of separating an image into two dif- 
ferent classes of pixels (generally the background and foreground, or dark 
and bright) by choosing the threshold that minimises the variance of the two 
classes. MODIFIED-CLEAN works best when small bright features are decon- 
volved before larger features, therefore successive rounds of Otsu thresholding 
were applied to the bright classes until the newest bright class contained a 
single pixel, or a maximum of 10 iterations. The single threshold that best 
selected small bright regions compared to the other candidate thresholds was 
chosen as t; via an “exclusivity”, €z, measure. 

When choosing t;, for each Otsu threshold, totsu, considered, €x is calculated 
using the fraction of selected pixels, fpix = Nbright/N, and the fraction of 
rejected range, frange = (totsu — min)/(max — min), where N is the number of 
pixels in the range of values, totsu is the computed Otsu threshold, Nprignt is the 
number of pixels above totsu, and min and max are the minimum and maximum 
values of the pixels in the range of values being considered. Combining these 
together gives the exclusivity as 


ez = (frange — foix)/(frange aF foix)- 


Of the candidate values of tots, considered, the one with the highest ey is 
chosen to be the initial threshold t; from which ts, is computed. 

Once the deconvolution is complete and the component delta-functions 
(CLEAN-components) are found, we do not re-convolve the final CLEAN- 
components with a “clean beam” as is standard practice. This is because one 
of the main goals of this deconvolution is to remove the “halo glow” effect 
of the adaptive optics, which significantly impacts limb-darkening Minnaert 
parameter estimations made from the original data. 

Before deconvolution, instrumental artefacts, missing data, and features 
smaller than the PSF were identified using an algorithm based on singular 
spectrum analysis (SSA) [39] and interpolated over. This reduced the number 
of iterations required and ameliorated one of the main problems with the family 
of CLEAN algorithms, their propensity to form speckles and ridges even with 
modifications [30]. 

The SSA artifact and scale detection algorithm was constructed to remove 
(if possible) or lessen (if not) the impact of instrumental artifacts and spurious 
noise on the MODIFIED-CLEAN algorithm. SSA is similar to Principle Compo- 
nent Analysis [40] in that it breaks down a dataset into m components, such 
that X = Xı +X2 +... +Xm = oy, Xi, for some dataset X, and components 
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are usually recorded in terms of decreasing relevance, ø (eigenvalue for PCA, 
singular value for SSA). However, SSA is fairly simple to extend to 2D data 
sets. Typically the SSA components of an image with large singular values cor- 
respond to the image’s subject, while a combination of components with small 
singular values corresponds to the noise in an image. The SSA of the PSF was 
used to determine the approximate lowest singular value, on, upon which some 
variation could be due to real data. Then for singular values lower than on, 
such that i > n, a cumulative density function F;(x) = P(x; < x) was con- 
structed for the values x;; in each component X;. Each pixel was then assigned 
a score derived from their deviation from the mean, S(x) = (2F;(x) — 1)?, 
and the scores for every pixel was averaged element-wise across the compo- 
nents S(j) = (1/(m—n)) 0", S;(a;;). Therefore, S(j) gives us a map of how 
consistently and how far the j*" pixel deviates from the median at each scale 
smaller than the scale of the PSF. Above a cutoff value, smax, a pixel was 
interpolated over. For this data we used Smax = 0.995. 

Deconvolution was continued until one of the following was true: 1) the 
RMS of the residual was reduced to 1% of the original image’s value; 2) the 
absolute brightest pixel was reduced to 1% of the original image’s value; 3) 
5000 iterations were performed. 

Lucy-Richardson (LR) and Maximum Entropy (ME) deconvolution algo- 
rithms were also explored as an option, but rejected. LR-deconvolution is 
very susceptible to “creating sources” /mottling in smooth regions of extended 
sources, and traditional methods of overcoming this (such as only applying 
LR above some radiance threshold or Tikhonov regularisation) did not align 
with our goal of correctly estimating the Minnaert limb-darkening parame- 
ters. ME-deconvolution methods, while widely used and considered accurate, 
are computationally expensive making their use across an entire datacube 
impractical. 

Other deconvolution techniques have been implemented to address 
the difficulties many algorithms have with extended sources. Differential 
deconvolution[41], for example, uses an estimated ‘background source’ image 
to remove extended regions, and thus change the problem into a more tradi- 
tional point-like source deconvolution. This technique and others like it [e.g., 
42] have definite advantages when working with extended sources that are pri- 
marily a flat(ish) background with bright point-like features. However, in our 
case they do not remove the extended source completely. The dark spot is itself 
an extended region of lower-brightness, thus even after accounting for the back- 
ground source we would still have an extended (albeit fainter) source to worry 
about which LR-deconvolution and the standard CLEAN algorithm would still 
struggle with. In the future, using differential deconvolution with MODIFIED- 
CLEAN as the base algorithm should be possible, and may be a fruitful avenue 
to pursue in further work. 


Data Availability. The raw VLT/MUSE datasets studied in this paper 
(under ESO/VLT program: 0104.C-0187) are available from the ESO Portal 
at https://archive.eso.org/eso/eso_archive_main.html. The reduced raw and 
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deconvolved ‘cubes’ for the observation IDs: 6 — 10, discussed in this paper, 
are available at https://doi.org/10.5281/zenodo.7594682. Data files associated 
with this analysis are available at https://doi.org/10.5281/zenodo.7620656 


Code Availability. The NEMESIS radiative transfer and retrieval code[17| 
used in this study is open-access and is available for download from GitHub or 
Zenodo[18]. The deconvolution code described in this paper is python-based, 
and still under development. However, the current version of this software is 
available from the corresponding author upon reasonable request. 
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Fig. 1 Processed Obs-6 VLT/MUSE observations. Panel A shows the observed and pro- 
cessed images at 551 nm, 831 nm, and 848 nm, with the first column showing the original 
MUSE observations, the second column showing the deconvolved and artefact-corrected 
images and the third column showing limb-darkening-corrected and contrast-enhanced 
deconvolved images. North is at top right. The application of our deconvolution algorithm 
gives a dramatic improvement in spatial resolution: in addition to NDS-2018, now just visi- 
ble at 551 nm at 15°N, we can see the South Polar Wave (SPW) near 60°S at 551 nm, and 
also a deep bright spot (DBS-2019) to the south west of NDS-2018, which is visible at 831 
nm, but not at 551 and 848 nm. Also visible in the 831 and 848 nm images are shallow bright 
spots (SBS) near the left-hand terminator, which are visible at all wavelengths longer than 
~ 600 nm. Panel B shows a false-colour image of the contrast-enhanced and limb-darkening- 
corrected data, where planetocentric latitudes and central meridian longitudes are spaced 
by 30°. Here the red channel is the 831-nm image, green is the 511-nm image and blue is the 
848-nm image. Panel B shows the close proximity between NDS-2018 (dark oval) and DBS- 
2019 (white/red spot to lower left of NDS-2018). Panel C compares the differences between 
the spectra observed in the NDS-2018, DBS-2019 and SBS regions (averaged over a 2 x 2- 
pixel box centred on the features) and the expected background spectra at these locations 
(See Methods), showing the very different reflectance signatures of these features. The errors 
on these spectra are set to the standard deviation of the values in the 2 x 2 pixel-averaging 
boxes and are shown for reference. N.B., the gap from 578 — 605 nm is set aside for the laser 
guide star. 
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Fig. 2 Analysis of NDS-2018 difference spectrum. Panel A shows the observed difference of 
the NDS-2018 spectrum from the ‘background’ spectrum expected at the NDS-2018 location 
from a Minnaert analysis of the 10—20°N latitude band, shaded with error limits set to the 
standard deviation of reflectivities in the 2 x 2-pixel averaging box. This difference spectrum 
is compared with two simulations: 1) modifying the opacity of the Aerosol-1 layer only 
(‘mod1’); and 2) modifying both the opacity and the imaginary refractive index (nimag) 
spectrum of the Aerosol-1 particles (‘mod2’). The gap in the spectra over 578 — 605 nm, 
indicated in grey, is set aside for the laser guide star. In these simulations we fitted to the 
observed NDS-2018 spectrum and also to the spectrum expected at the same location from 
a limb-darkening analysis of the observations in the 10 — 20°N latitude band, and show 
the difference in the fitted reflectivities. As can be seen, modifying only the opacity of the 
Aerosol-1 layer (‘mod1’) provides a poor match with the observed difference spectrum, while 
modifying both the opacity and scattering properties (‘mod2’) leads to an excellent fit. The 
fitted and assumed njmag spectra (shaded to indicate the retrieved error limits) are shown 
in Panel B. Here the nimag spectra of the Aerosol-2 and Aerosol-3 particles are derived 
from the analysis of the disc-averaged spectra, while the njmag spectrum of the Aerosol- 
1 particles is re-retrieved from NDS-2018 and background spectra. The vertical profiles of 
aerosol opacity (opacity/km at 800 nm), again shaded to indicate the formal retrieval errors, 
are shown in Panel C. The vertical structure of the Aerosol-2 and Aerosol-3 particles is 
determined from a Minnaert analysis of the 10 — 20°N latitude band, while the Aerosol-1 
opacity is re-retrieved from the NDS-2018 and ‘background’ spectra. In Panel B, it can be 
seen that in the dark spot nimag has increased for Aerosol-1 for A < 700 nm, lowering the 
single-scattering albedo. 
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Fig. 3 As Fig. 2, but modelling the observed DBS-2019 difference perturbation spectrum. 
Panel A shows the observed difference of the DBS-2019 spectrum from the background, 
compared with a simulation where we have either: a) modified the opacity and the nimag 
spectrum of the 0.1-m mean radius Aerosol-1 particles (‘mod3’); or b) increased the mean 
radius of the Aerosol-1 particles from 0.1 to 0.7 um in the spot and also modified the opacity 
and Nimag spectrum (‘mod4’). Adjusting the Aerosol-1 particles in either way leads to good 
fit at all longer wavelengths, but a slightly better fit at shorter wavelengths is obtained when 
we also increase the Aerosol-1 mean particle radius. The assumed and perturbed aerosol 
Nimag Spectra and vertical profiles of aerosol opacity (opacity/km at 800 nm), shaded to 
include the formal retrieval errors, are shown in Panels B and C, respectively. 


